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Simulation methods based on stochastic realizations of state vector evolutions are commonly 
used tools to solve open quantum system dynamics, both in the Markovian and non-Markovian 
regime. Here, we address the question of waiting time distribution (WTD) of quantum jumps for 
non-Markovian systems. We generalize Markovian quantum trajectory methods in the sense of 
deriving an exact analytical WTD for non-Markovian quantum dynamics and show explicitly how 
to construct this distribution for certain commonly used quantum optical systems. 



I. INTRODUCTION 

An open quantum system interacting with its envi- 
ronment undergoes non-unitary evolution and typically 
loses its quantum properties, such as entanglement, due 
to decoherence P] . Whilst the theory of Markovian dy- 
namics in terms of semigroups and completely positive 
trace preserving maps is fairly well understood since the 
pioneering work of Lindblad, Gorini, Kossakowski, and 
Sudarshan [51 [5] , non-Markovian quantum dynamics dis- 
playing memory effects has become under active study 
during the recent years. The advances here include the 
development of simulation schemes [IHTB] , the limits for 
the existence of physically valid dynamical maps |17j . 
the discussion about the applicability of different types 
of master equations [TB], the very definition and quan- 
tification of quantum non-Markovianity |19H21) . and the 
role of initial correlations between the system and its en- 
vironment p5H?7] . Moreover, it is also possible to control 
and quantify experimentally the non-Markovian features 
of quantum dynamics 28, 29J and the influence of initial 
system-environment correlations [30l [31] . Subsequently, 
this progress allows to look for ways how non-Markovian 
features with memory effects can be exploited for quan- 
tum information processing |32) , and for quantum control 
and engineering tasks [33l |34| • 

Here, our focus is on fundamental aspects of non- 
Markovianity and, in particular, on the jumplike stochas- 
tic unravellings, or simulation schemes, for open system 
dynamics [5l [9Hl4| , [35 | [36] . For Markovian systems, some 
of the most popular stochastic schemes include the Monte 
Carlo Wave Function (MCWF) [371 131] and Quantum 
Trajectory (QT) [39H4T] methods. In both of these meth- 
ods the time evolution of a singe realization consists of 
periods of continuous deterministic evolution interrupted 
by stochastic jumps, i.e., both methods simulate a piece- 
wise deterministic stochastic process (PDF). In MCWF 
method, the time evolution of a single realization pro- 
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gresses in a stepwise fashion, e.g. during each time step 
we decide whether the realization evolves deterministi- 
cally or jumps. The mean time evolution of the ensem- 
ble of realizations, over small time increments, matches 
with the solution of a Markovian master equation for the 
density matrix (for the first order in time increment). 
The central concept for the QT methods, in turn, is 
the waiting time distribution (WTD). The random jump 
time of the realization can be sampled from the WTD 
and the state vector is directly evolved deterministically 
till this point. Solution to the Markovian master equa- 
tion is formed from the weighted average over all possible 
stochastic evolutions that realizations might take. Gen- 
erally speaking, the MCWF method exploits the incre- 
ments of the WTD while the QT uses the full exact form 
of the WTD. 

A few years ago MCWF was generalized to non- 
Markovian region by Non-Markovian Quantum Jump 
method (NMQJ) [111 [SI [35]. In NMQJ, evolution of 
the ensemble average over a time step St, matches with 
the solution given by local in time master equation with 
possibly temporarily negative rates. The central ingredi- 
ent of the NMQJ method is a quantum jump which can 
restore coherence, e.g., by returning the stochastic real- 
ization to the superposition which was destroyed earlier. 
Formally, the probability of the reverse jump can be cal- 
culated using the concept of positive definite jump prob- 
ability density |36| . However, to the best of our knowl- 
edge, the QT methods - without using the auxiliary ex- 
tensions of the state space of an open quantum system - 
have not yet been extended to the non-Markovian region. 
The main obstacle here has been the fact that the WTD 
for non-Markovian systems, when calculated along the 
Markovian line of reasoning, displays oscillations which 
render its physical meaning invalid and prevent the tech- 
nical implementation of the simulations, whilst the math- 
ematical calculation of the WTD still is, in some sense, 
correct. 

With the help of the insight provided by the NMQJ 
method and the concept of positive definite jump proba- 
bility density, we derive a general analytical form of the 
waiting distribution, which is both physically and mathe- 
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matically correct for non-Markovian quantum dynamics. 
This is the main result of our paper. We thereby general- 
ize the QT formalism into the non-Markovian regime and 
show explicitly how to construct the WTD for some com- 
monly used quantum optical systems. It is worth keep- 
ing in mind here that, as already featured in the NMQJ 
method, the stochastic realizations depend on each other 
as a consequence of the memory effects. Moreover, it 
has been recently argued that non-Markovian unravel- 
ings can not be interpreted as stemming from continuous 
measurement of the environment |42j despite of some at- 
tempts in that direction [331 H?]. It seems to us, that 
the functional form of the derived non-Markovian WTD 
indicates the former choice of answers. 

The structure of the paper is the following. In Sec. [TTj 
we introduce the PDP corresponding to the NMQJ 
method and most importantly the positive definite jump 
probability density. In Sec. |III[ we give the general form 
of the waiting time distribution and connect it to the 



PDP defined in Sec.|ll] In Sec.|TV|we present some quan- 
tum optical examples that illustrate the general construc- 
tion of the WTD and the effects of non-Markovianity, in 
Sec. [Vjwe present some further discussion about our re- 
sults and conclude in Sec. IVII 



II. 



PIECEWISE DETERMINISTIC PROCESS 
FOR NON-MARKOVIAN SYSTEM 



In this section we formulate non-Markovian piecewise 
deterministic process for pure states tp [131 El ISS] . Re- 
duced state of the system, p, is obtained as an ensemble 
average 



p(t)=E[|V')(V'|]= / d^F[^,t]|^)(^|, 



(1) 

where d^ = DipD-tjj* is a singular volume element of 
the Hilbert space of the system and P[ip,t] is time de- 
pendent, phase invariant probability density functional 
concentrated on the surface of a unit sphere \ — 1). 
p{t) solves also the following time convolutionless (TCL) 
master equation [T] 

+ X{t) (c.p{t)Cj - ^ {p{t), cjc,}^ 
^~ih-^[Hs{t),p{t)] 

+ 5]A+W (c,p{t)Cj \ [p{t\C]C,]^ 
-Y^A^it) (c,pit)Cl - ^ {pit), Clc}') , (2) 



where Ai{t) is time dependent decay rate. After the sec- 
ond equality sign, we have split the decay rates into two 
components A^(t) = (|Aj(i)| ± Aj{t)) /2 to account bet- 
ter for the overall sign of the decay rate [TJ. However, 



note that A~^{t) are non-negative for all times t. From 
now on we assume that h = 1. Operators Cj are called 
jump operators and we make a simplifying assumption 
that they are time invariant. Unnormalized states are la- 
beled with ip and normalized with ip. We formulate the 
process for pure initial states only, since mixedness adds 
no novelty here. 

Between two subsequent jumps at times T and t ^ T + 
T (t > 0), pure states evolve deterministically according 
to an effective non-Hermitian Hamiltonian 



such that state tp(t) is expressed as 
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(3) 



(4) 



where V't(''") satisfies the Schrodinger equation 4't{t) = 
~iHcs{T -\- t)^t{t) with the initial condition '(/'t(O) = 

Discontinuous part of the process consists of jumps be- 
tween different pure states. Given that the process is in 
pure state ^p, conditional jump probability density from 
a source state ip to a target state (p using a channel k 
during a time interval [T, T -I- St] is [35] 



Pkm.T] ^6tA+{T)\\Ck^{T)f5 
P[<P,T] 



+ 5tAl{T) 
X 5 (i^{T) 
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Above, (5-functional satisfies J d(p6{ip — (p)F[(p] = F[ip] 
where F is an arbitrary smooth functional. The 6- 
functionals in Eq. ^ give temporal channel-wise stochas- 
tic connection between different regions of projective 
Hilbert space (global phase of the states is irrelevant). 
Connection of the positive part (i.e. part proportional 
to A^) is of one-to-one type: Tp which corre- 

sponds to Markovian quantum jumps. Interestingly, con- 
nection of the negative part is one-to-many type: Each 
source state "0 may jump to one of the states {cp} that 



satisfy ^ 



Ck4> 

IIC.0II 



provided that the corresponding jump 



probability is nonzero. It follows that the connection 
provided by the negative part requires the knowledge of 
the different states in the pure state decomposition of p, 
since the range of the one-to-many mapping is not ob- 
tainable from the structure of Eq. To summarize, a 
negative channel induces a one-to-many mapping for the 
pure states and therefore, in general, one decay channel 
connects several different regions of the projective Hilbert 
space stochastically. 

Next we sketch the stepwise progression of the PDP, 
more details may be found in Refs. [351 136j . During an 
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interval / = [T+T,T + T + 5t], a realization of the process 
in state ip may either jump or evolve deterministically. 
The total jump rate away from state ^ during the interval 
/ is the total jump probability to any other state via any 
channel divided by the length of the interval 



r[^,T- 



1 

Tt 



(6) 



Therefore, with probability 1 — T [ip,T + t] 6t, the real- 
ization does not jump away from state ijj but evolves 
deterministically. Deterministic evolution is governed by 
the Schrodinger equation and Eq. ([s]). With probabil- 
ity T [ip^T + t] St the realization jumps; the target state 
of the jump is chosen from the probability distribution 
^ri^jj^T+T^st ■ After the stochastic evolution of the ensem- 
ble over a small time step, the average over the ensemble 
provides us Eq. ^ for the first order in St. 



III. WAITING TIME DISTRIBUTION FOR 
NON-MARKOVIAN SYSTEM 

In this section we derive the general form of the waiting 
time distribution, which is valid also for non-Markovian 
systems, starting from the positive definite jump proba- 
bility density. We also provide a formula for estimating 
the WTD from a sample of realizations. 



A. Analytical WTD 

By definition, the waiting time distribution F{t\'iJj,T) 
is a conditional probability distribution function which 
gives the probability for the next jump to occur during 
a time interval [T, T + t] conditioned on that at time T 
the state of the realization is known to be ^/^ [T]. 

Probability for a jump to occur during a short time 
interval / — [T + t,T + t + St] away from state ip is 
then SF{T\i;,T) = F{t + St\'il;,T) - F(t|V',T), which is 
equal to the probability of having no jumps before T + t 
and a jump during the following St, i.e. SFlrlip^T) — 
(1 - F{T\'i(j,T))T[i},T + T]St. Dividing both sides by St 
and taking the limit St we obtain the following dif- 
ferential equation that every valid WTD must satisfy [T] 



dr 



F{T\i^, T) = [1 - F{t\^, T)] r[^,T + T]. (7) 



This can be solved formally with an initial condition 
F(0|V',T) = 0, such that 



i^(r|i^,T)=l-exp|-^ ^ clsr[^,s]| 



(8) 



Then, by using Eqs. ([s]), ([6]), and Eq. ([s]) we can obtain 
the following form for the generic WTD corresponding to 



Eq. ([2|: 

F(r|V',r) =1 -cxpj ds J dq 
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Terms proportional to depend only on the state of 
the particular realization, its deterministic time evolution 
and the quantities obtainable from Eq. Terms pro- 
portional to are more complicated, since they depend 
on the probability functionals and on the deterministic 
time evolution of other states to which the realization 
might jump via a channel- wise one-to-many mapping. 

Random waiting time t* is sampled from waiting time 
distribution by comparing a random number rj to the 
WTD: T*{r]) = min{T|F(r|7/>, T) > r]} [J. Probabilities 
P[ijj,s] appear on the right hand side of Eq. ^ and they 
are modified each time a jump occurs in the ensemble. 

When all decay rates Ai{t) for all times t are non- 
negative in Eq. (pi), the total jump rate away from pure 



state ?A is T[ip,t\^ J2k ^kit)\\Ckipit)\\'^ . Inserting this 
into Eq. ^ and taking into account the deterministic 
evolution of ■0, we obtain the following familiar Marko- 
vian limit for the WTD [2 [Ml iOl SS] : 



F{t\^P,T) 



WMOW 



(10) 



Details of the derivation of the Markovian limit can be 
found in the Appendix [X] 

B. Estimation of WTD 

We assume that the reduced state of a non-Markovian 
open quantum system can be expressed at all times, as 
a linear combination of a finite number of, in general, 
non-orthogonal pure state projectors. Then we can write 
Eq. ^ as 



(11) 



Assume that we have a sample of Ns realizations from 
PDP in Sec. [TT| over a time interval [^oi^s] divided into 
Nt time steps. The samples are collected to an Nt x TVs 



matrix M where the element Mi 



means that a 



realization j is in state 4'^ at time ti = {i — l)St + to. 
Set of column indices /f of row i of M, give the indices 
of the realizations which are in state /3 at time ti. Hence 
each set has Ns elements, where the fcth element is 1 
if realization k is in state B at time ti, otherwise the fcth 
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FIG. 1. (Color online) Schematic figures of the different sys- 
tems, a) Two level system, b) A system and c) Ladder system. 



element is 0. |/f | = Y^^li{-^i)k is the total number of 
realizations in state /3 at time i. 

If we know that the realization r is in state a at time i^, 
then the discrete sample estimate for the probability to 
jump away from state a during the discrete time interval 



Wr{tkKa)^l- ^ 



-1 ' ' 



l=i+l 



(12) 



Naturally we have that Wr{ti\ti, a) = 0. The meaning of 
this equation is that the intersection of two sets consists 
of the indices of those realizations that were in a state 
a at the previous time and are still there at the present 
time. The number of such realizations is divided by the 
number of the realizations in a at time ti (beginning of 
the time interval). 



IV. CONSTRUCTION OF WTD FOR 
QUANTUM OPTICAL SYSTEMS 

In this section we construct the waiting time distri- 
bution explicitly for a few simple quantum optical sys- 
tems interacting with a leaky cavity mode. In Fig. [TJwe 
have presented schematically the different systems that 
we shall study. 



A. Two level system 

Orthonormal basis for the Hilbcrt space of the system 
is {|0), |1)}, where |1) is the groimd state and |0) is the 
excited state of a two level atom (TLA), (see Fig. [T^)). 
The initial state is i/'°(io) = co(io)|0) + ci(to)|l), which 
is the only state with non-trivial deterministic evolution. 
The state of the system is decomposed for all times t 
as p[t) = Po{t)\il^°){'4j'^\+ Pi{t)\ip^){il)\ where tp^ is the 
ground state. The detailed description of the system is 
given in Appendix |B 1[ 

The non-Hermitian Hamiltonian generating the deter- 
ministic pieces of the time evolution is obtained from 



Eq. (§ (see details in the Appendix IB 1|) . The total rate 
away from the deterministic state ilrjt) \s 



0, 



A(i) > 0, 
A(i) < 0, 



and the total rate away from the state -0^ is 



Jo, AW>0, 
\|A(i)|§g||CV'°(i)lP, AW<0. 



(13) 



(14) 



Inserting the rates ( 13 1 and ( |14[ ) as well as the analytical 
solutions of Appendix B 1 for the probabilities Pq (i) and 



Pi{t) into Eq. ([7f, we may solve a formal expression for 
the WTD. The solution depends on the particular path 
that one realization might take. For example, WTD is 
different for a jump ■0'^ — >■ -0^ somewhere in the inter- 
val [T, T + t] if the realization has made zero or two 
transitions before time T. We illustrate this in Fig. [5] 
where we have plotted the decay rate A(t), three sam- 
ple realizations and the WTDs for each realization solved 
from Eq. (It]) and also from Eq. (12). The initial state is 
'0°(O) = |0) an d we use parameter values 70 = 5A, (5 = 8A 
(see Appendix |B l\ and a sample size of 10^. Points of 
discontinuity in the waiting time distribution in panel 
e) of Fig. [2] correspond to jumps and since the state of 
the realization changes, the waiting time distribution also 
changes. We see that during periods of negative decay 
rate, the derivative of WTD is zero for realizations that 
are in state 0", since the jump rate is zero. 

This system is the simplest one since it has only one de- 
cay channel and the pure state decomposition of Eq. ( 11 ) 



consists of two states. Jump paths between the different 
states in the pure state decomposition show that, both 
the positive and the negative channel act as a one-to-one 
map in the projective Hilbert space of the system. 

It is interesting to consider the WTD for a realization, 
which jumps at some time during the first positive decay 
rate region and then makes a reverse jump during the 
first negative region. For the first positive region [ioi^i) 
we obtain 



J^(r|0°,to) 



.ll<(0)|p-||<(r)|| = 



(15) 



and for the first negative region [ti,t2) 



F(r|^i,ti) 



(16) 



When comparing the WTD of Eq. ( 15 ) for the positive 



jumps to the WTD of Eq. ( 16 1 for the negative jumps. 



we see that they are complementary: in the numerator 
the norm decrease of the state ■0° in the positive region 
is switched to a norm increase in the negative region and, 
the denominator in the negative region is the complement 
of the denominator in the positive region. 
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FIG. 2. (Color online) Initial state is 1 1/'°(0)) = |0) and param- 
eters are 70 = 5A, 5 = 8A and sample size was 10^. In panel 
a) we have the decay rate, in 6) — d) we have three different 
realizations. In panel e) we have WTDs for the realizations. 
Line style and color coding match with the sample realization. 
Lines are for exact numerical solution and markers for sample 
estimate. 



Equations ( 15 ) and ( 16 1 provide a simple way of doing 
a simulation for the TLA. For example, during the A:th 
negative region we could calculate the random waiting 
time using Eq. (16 1 with substitution t\ — t^, for each 
realization that are in the ground state. During the /cth 
negative period realizations that are not in the ground 
state do not have a possibility to jump. During the fcth 
positive period we would use Eq. ( 15 1 with ~^ ifc-i for 
jumps away from the state 'tlP . 



B. A-system 



Eq. ([3]) and Appendix B2). The state of the sys- 
tem p(t) can be decomposed for all times t as p{t) = 
Yi=^Pk{t)\'^^){^% states V'' = with /c = 1,2, are 
the ground states of the system. The probabilities ap- 
pearing in the decomposition are explicitly calculated in 
the Appendix [B2| 

Since we have two decay rates, we have four possible 
combinations of the decay rate signs. For each pure state 
of the decomposition we only present the decay rate sign 
combinations which lead to a non-zero jump rate away 
from the state under consideration. Other sign combina- 
tions would produce zero rate. Jump rate away from the 
state V'° is 



A,(i)||Q^"(i)|P, 



Ai, A2 > 0, 
A; > A A,- 



< 0. 



Jump rate away from the ground state 1/''^ is 

^[v;^^]=|A,(^)|^||CfcV°WII^ 



(17) 



(18) 



when Afe(t) < 0. This is the case irrespective of the sign 
of the other decay rate. 

Both channels, irrespective of the sign of the decay 
rate, are one-to-one maps. However, when both chan- 
nels are positive, i^P may be mapped to i^)^ or ■0^ when 
considering the effect of both channels. All rates are pro- 
portional to ||Cfe?/'"WlP = IcoWP- 

Some realizations are plotted together with their WTD 
in Fig. [3] The initial state is '(/'^(O) = |0) and we use 
parameter val ues 7 n"^'^^ = 5A, (5^^^ = 4A, (5^^^ — 8A 
(see Appendix B2) and a sample size of 10^. We ob- 
tain an interesting expression for the WTD for a reverse 
jump 7/;^ — >■ ilP , if we let Ai(t) < during time inter- 
vals [sJjSj], [■53,54], etc. If a jump to the state i])^ oc- 
curred at time T e [^o, s}] then the probability for jump 
away from i])^ somewhere in the interval [to, T-l-r], where 



[4n- 



F(r|V'i,T)-l 



Pi(si)Pi(4)''>i(sL_i)- 



(19) 



Since Ai(t) < when t S [s2n-i' ■52n]> probabilities 
P^^^\n) < Pi{^2n-i)- Therefore, each fraction is smaller 
than unity and F{t\iP^,T) is monotonically increasing 
function. 



Let us indicate the basis for the Hilbert space of 
the system with {|0), |1), |2)}, where |1) and |2) are 
the ground states and |0) is the common excited state. 
Schematic representation of this system is in Fig. lip). 
Initial state is V°(<o) = co(to)|0) + ci(<o)|l) + C2{to)\2), 
which is the only state with non-trivial deterministic evo- 
lution (see Appendix B 2 for more details). 



C. Ladder system 

We label the orthonormal basis for the Hilbert space 
of the system with {|0), |1), |2)}, where |0) is the excited 
state, |1) is the middle state and |2) is the ground state. 
Schematic representation of this system is in Fi, 



Deterministic evolution is generated by Hcs{t) (see The initial state is of the form '>p (to) = co{t, 
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FIG. 3. (Color online) Initial state is |iA°(0)) = |0) and pa- 
rameters are 7g^'^' = 5A, ij'^' = 4A, (5'^-* = 8A and sample size 
was 10^. In panel a) we have the decay rate, in 6) — d) we 
have three different realizations. In panel e) we have WTDs 
for the realizations. Line style and color coding match with 
the sample realization. Lines are for exact numerical solution 
and markers for sample estimate. 



ci(io)|l) + C2(to)|2)- The deterministic evolution is gen- 
erated by Hce{t) (see Appendix |B 3| and Eq. ([3|). For 
all times t the state of the system p{t) may be decom- 

posed as pit) = ELo^fcWlV'')(V'1, where = 
with k — 1,2, are the middle and the ground states, 
respectively. Analytical expressions for the probabilities 
Pi (t) are in Appendix B 3 For this system, the o nly s tate 
invariant in respect to H^e is "0^ (see Appendix B3). 



As in Sec. |IVB| we write down only those combinations 
of the decay rates that give a non-zero jump rate. For 
the initial state tp'^{t) we have 



and for the middle state ip^{t) we have 



^A2(t), 

A2W 

+|Ai(t)|^||CiV°WlP, 
lAiWI^IICiV-^WlP, 



Ai,A2 > 0, 

A2 > A Ai < 0, 
Ai,A2 < 0, 

(21) 



and for the ground state ip'^ we have 

^[v.^^]=|A2(^)|(^||C720°(^)ll^ 



Pijt) 

P2{t) 



(22) 



when A2{t) < irrespective of the sign of Ai(t). 

Channel 1 maps ■0° to -0^ and channel 2 maps ■0" to -0^ 
and to ip"^. When decay rates are negative, channel 
1 maps "0^ to 0*^. However, channel 2 maps "0^ to ip'^ or 
when negative. Therefore, when jump to channel 2 
occurs when it is negative, we still have a probability dis- 
tribution over the two different target states from which, 
we have to choose the actual target state for the jump. 

In Fig. [i] we have plotted decay rates and three realiza- 
tions with their respective WTD. There, the initial state 
we use is i/'°(0) — |0) and the param eters are 79^'^'' = 5A, 
<5(i) = 8A, (5(2) = 4A (see Appendix B3) and we used a 
sample size of 10^. 

It has been shown in Ref. |35j that for some parameter 
values the approximations made while obtaining the mas- 
ter equation for this level geometry fail, which is mani- 
fested by the breakdown of positivity. This is due to 
the fact that the population of the ground state 0^ is 
drained completely while the decay rate A2{t) is still 
negative. This causes Eq. ( |22| ) to diverge. Let us as- 
sume that A2{t) < during intervals the [^3,^4], 
etc. and that at time T € [^oi^i) the realization jumps 
to state ■0^. Assuming that T + t € [^in-i'^ln] then, the 
analytical form for WTD reads 



F{r\r 



' ' P2{tl)P2{tl) 



P2{T + 
P2{tL- 



(23) 



From Eq. (23) we see that if limt/_j.j'+T- P2{t') — 0, wait- 



ing time distribution reaches unity in finite time but it 
is still well defined. Dynamical consequences of this are 
that a simulation method utilizing full WTD would not 
break down. Instead, population of the state -0^ would go 
to zero and the total population is distributed between 
the pure states ip^ and . 



V. DISCUSSION 



E,A,(OI|Cfc^°WIP, 

A,(t)||av"(t)|P, 



Ai,A2 > 0, 
Ai > A Aj < 0, 
(20) 



The positive definite jump probability density of 
Eq. ([5| shows that there is a correlation between the 
different regions of the projective Hilbert space. There- 
fore, a general form of the WTD in Eq. ^ is complicated 
since it takes the correlation into account cumulatively. 
On the other hand, it confirms that the realizations of 
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considered in this work this happens in the TLA always; 
in the A-system always for ground states and for the state 
-0'^ , when the decay rates have the same signs; and in the 
Ladder system for the ground state always, for the middle 
state when the decay rates have the opposite signs and 
for the state "0" when the decay rates rates have equal 
signs. 

For a short time interval 5t we can approximate the 
full WTD as 



F{5t\'^,T)^V[^{T),T]5t 



mv (25) 



Thus, for a short time interval the total jump rate is 
resolvable in (channel, target state) pairs: each channel 
maps a source state to a target state (one-to-one rela- 
tion for the Markovian jumps and one-to-many for the 
non-Markovian jumps). During this short interval, the 
occurrence of a jump excludes the possibility of another 
jump at the same interval to another channel. In WTD- 
based methods these individual contributions are cumu- 
latively gathered together. The process may be reset 
after any time interval At, after which a new random 
number must be drawn. In the limit At — > (5t, stepwise 
method emerges. 



FIG. 4. (Color online) Initial state is |V'"(0)) = |0) and pa- 
rameters are 7g^'^' = 5A, ij'^' = 8A, (5'^-* — 4A and sample size 
was 10®. In panel a) we have the decay rate, in 6) — d) we 
have three different realizations. In panel e) we have WTDs 
for the realizations. Line style and color coding match with 
the sample realizations. Lines are for exact numerical solution 
and markers for sample estimate. 

the PDF considered in this paper do not form a trajec- 
tory, i.e. continuous measurement interpretation can not 
be necessarily made. This happens because it is not pos- 
sible to express the WTD for a given realization in terms 
of that particular realization only. This is the argument 
used already by Gambetta and Wiseman in the context 
of non-Markovian quantum state diffusion [5] but it can 
be also applied here. For further discussion on this highly 
non-trivial topic, we refer the reader to Refs. [i^ El ? ]. 

In the case that there is a state 4'^ in the pure state 
decomposition of p{t) that acts only as a source state 
for jumps for some period [T, T -f r], then the WTD is 
quite simple over this period. During this period, the 
probability of the state ip^ in the pure state decom- 
position changes only by jumps away from that state. 
Hence, we have the following identity Pfc(T + t) ~ 
Pk{T) ~ F{t\\1}^ ,T)Pk[T) from which we can solve 

n.|*'.r).^MI±I>^, (24) 

where Pk{T) 7^ is assumed. In the examples that we 



VI. CONCLUSIONS 

We have derived a general waiting time distribution 
of quantum jumps for open quantum systems following 
non-Markovian dynamics. In this sense, our results gen- 
eralize the QT methods into the non-Markovian regime. 
The distribution is a well defined conditional probabil- 
ity distribution function which takes into account in a 
proper manner the bidirectional probability flow between 
different regions of the projective Hilbcrt space of the sys- 
tem. The WTD includes probabilities which are present 
in the pure state decomposition of the reduced system 
state, i.e., the realizations of the process depend on each 
other - a feature stemming from the memory effects and 
present already in the NMQJ method. Our results seem 
to confirm the view that the realizations of the PDF, that 
the WTD govern, do not form a trajectory, therefore the 
FDF can not be interpreted in terms of a continuous 
measurement of the environment. We have constructed 
the WTD explicitly for some quantum optical systems 
and also discussed the cases when the calculation of the 
WTD can be simplified. 

Our work complements the theory of Monte Carlo 
methods for non-Markovian systems and the WTD con- 
cept familiar from Markovian regime is now also well- 
defined for non-Markovian systems. We hope that this 
work stimulates further research for non-Markovian dy- 
namics and especially inspires new directions in the de- 
velopment of simulation tools for open quantum systems. 
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Appendix A: Markovian limit 

In the Markovian limit, decay rates Ai{t) > 0, for all 
times t. Then, the jump rate away from state ^j) at time 
T + s is 

r[ij,T + s] =J2\{T + s)\\CMT + s)\\' 

i 



The time convolutionless master equations for the ex- 



his) 



(Al) 



where we have used notation of Eq. Q. On the other 
hand, 

^J\Ms)\\' = -J2^^('^ + ')\\^^Ms)\\', (A2) 

i 

where we have used Schrodinger equation and Eq. (|3]). 
Therefore, the total jump rate is 



r[^,r + ,s] = -A^^^!^^-Ai,||^^(,)l|2. (A3) 

Now, using Eq. (Is]) we obtain 



F(T|V,T)=l-exp<j 1^ ds\n\\Ms) 



-1 — exp < In ■ 



^ ^IIV-tWIPJ 

^ ||^T(0)|p-||^T(r)|p 

IIV^t(0)|| 

Appendix B: System definitions 



(A4) 



ample systems in Sec. IV are all special cases from the 
following general form 



-f ^Afe(t) (^Ckp{t)Cl~^ {pit),ClCk}^ , (B2) 

where Sfe(i) is the time dependent Lamb shift, Afe(t) are 
the time dependent decay rates and Ck are the time in- 
dependent jump operators. For simplicity, we have as- 
sumed in the actual calculations that Sk{t) — 0. We use 
TCL4 approximated analytical form for the decay rate 
[1] corresponding to a spectral density of Eq. (Bl) 



A{t) 



1 - e 



7oA^ 

6 



-At 



cos((5t) — — sin((5t) 



1 - 3 



(e^* - e"^* cos(2,5t)) 



3- 



1 - ( - ) Xtcos{6t)+4: 1 + i^jj dtsm{St) 
e-^*sin(2(5i)|. (B3) 

1. Two level system 



Master equation: 

p(t) = -z.(i)[|0)(0|,p(t)]+A(t)|l)(0|p(i)|0)(l| 
1 



-A(i){pW,|0)(0|} 



(B4) 



Jump operator: 



C=.|1)(0|. 



(B5) 



Solution for the probabilities in the pure state decompo- 
sition of Sec. IIV AI 



Poit)=\\i'lit~toW, 
Pi{t) =l-||<(i-to 



(B6) 
(B7) 



We are considering two- and three level atoms inter- 
acting with a leaky cavity mode. The spectral density of 
the cavity is 



1 



7oA2 



27r (w - Wc)2 + ^2 ' 



(Bl) 



where 70 is the coupling constant, A is the width of the 
Lorentzian and Wc is the cavity resonance frequency. An- 
other important parameter is the detuning of the atom 
from the cavity resonance: 5 — uja — lOc, where uja is one 
of the transition frequencies of the atom. 



2. A-system 

Master equation: 
pit) - - zsi(t)[|0)(0|,p(i)] - is,{t)[\0){0\,p{t)] 



Ai(t) 

A2{t) 



|l)(0|p(t)|0)(l|--{p(t),|0)(0|} 
|2)(0|p(0|0)(2|-i{p(t),|0)(0|} 



(B8) 
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Jump operators: 



3. Ladder system 



Master equation: 

p{t) = - iSi(t)[|0)(0|,p(i)] - *S2(i)[|l)(l|,p(t)] 



C^i=|l>(0|, 
C2=|2)(0|. 



(B9) 
(BIO) 



+ Ai(t) 

+ A2(t) 

Jump operators: 



|i)(o|pW|o)(i|-^Mt),|o)(o|} 

\2){l\p{t)\l){2\^\{p{t),\l){l\} 



Solution for the probabilities in the pure state decompo- 
sition of Sec. ITV^ 



Gi = \l){% 
C2 = |2)(l|. 



(B12) 



(B13) 
(B14) 



P,{t)= / dsA,(s)|5o(s)|^ 
Jt„ 



(Bll) 



Solution for the probabilities in pure state decomposition 
of Sec. HVCl 

P,{t)=Ul{t-t,)\\\ (B15) 

P,{t) ==e-^=W rdsAi(.)e-^^(^H^^(^)||CiV^to)|^ 

(B16) 

P2(i) = (l-e~^^(*^) ||C2^"(to)||'+ /' dsA2(s)P2(s), 

(B17) 



to 



for J e {1,2}. 



where Di{t) — ds Ai(s). 
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